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Abstract 

We present a lattice Boltzmann solution of the equations of motion de- 
scribing the spreading of droplets on topologically patterned substrates. We 
apply it to model superhydrophobic behaviour on surfaces covered by an ar- 
ray of micron-scale posts. We find that the patterning results in a substantial 
increase in contact angle, from 110° to 156°. The dynamics of the transition 
from drops suspended on top of the posts to drops collapsed in the grooves 
is described. 

1 Introduction 

From microfluidic technology to detergent design and ink-jet printing it is vital to 
understand the way in which drops move across surfaces. The dynamics of the 
drops will be affected by any chemical or topological heterogeneities on the sur- 
face. Until recently such disorder was usually regarded as undesirable. However 
with the advent of microfabrication techniques it has become possible to control 
the chemical or topographical patterning of a substrate on micron length scales, 
leading to the possibility of exploring new physics and to novel applications. 

A beautiful example of this, inspired by the leaves of the lotus plant, are su- 
perhydrophobic substrates. These are surfaces which are covered with posts on 
length scales of order microns. As a result of the topological patterning they are 
strongly repellent to water droplets which show contact angles up to 160° JHElEll. 
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This should be compared to more traditional ways of increasing the contact angle, 
surface coatings and chemical modifications of the substrate, where it is difficult 
to achieve an angle of more than 120°. Superhydrophobic substrates have many 
potential applications, for example as materials for raincoats or windscreens. The 
evolutionary advantage to the lotus appears to be the easy run-off which helps to 
clean the leaves of the plant. 

Droplets on a superhydrophobic surface can be in two states, suspended where 
the drop sits on top of the posts with pockets of air beneath it, or collapsed where 
it wets the grooves between the posts. Several authors have calculated the free 
energies of the suspended and collapsed states using approximations based on the 
Cassie-Baxter and Wenzel laws respectively 00. They have shown that both 
states can be fhermodynamically stable with the phase boundary between them 
depending on the surface tension and substrate geometry. It has also been argued 
that the suspended drop is often observed as a metastable state as it has to cross a 
free energy barrier to fill the grooves. However the kinetics of the transition to the 
collapsed phase is not understood: it is not accessible to the equilibrium theories 
and it is hard to probe experimentally. 

Therefore in this paper we present a numerical solution to equations which 
are able to describe both the static and dynamic behaviour of droplets on topo- 
logically patterned substrates. The droplet dynamics is described by the Navier- 
Stokes equations for a liquid-gas system. Its equilibrium behaviour corresponds to 
a chosen free energy functional so that appropriate thermodynamic information, 
such as the surface tension and the contact angle, are included in the model. We 
solve the equations of motion using a lattice Boltzmann algorithm. This approach 
has a natural length scale, for fluids such as water, of order microns where much 
of the exciting new physics is expected to appear. The method has already shown 
its capability in dealing with spreading on surfaces with chemical patterning 0. 

In section 2 we summarise the algorithm and, particularly, describe the new 
thermodynamic and velocity boundary conditions needed to treat surfaces with 
topological patterning. In section 3 we present results for a substrate patterned by 
an array of posts and show that the patterning leads to a substantial increase in 
contact angle in agreement with experiments. We then explore the kinetics of the 
transition between the suspended and collapsed droplet states. Finally we discuss 
directions for future work using this approach. 
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2 The mesoscopic model 



2.1 Equilibrium free energy 

We consider a liquid-gas system of density n(r) and volume V. The surface of the 
substrate is denoted by S. The equilibrium properties of the drop are described by 
the free energy 

* = Iv (^ (n) + ^ dV + Is 0) 

where Einstein notation is understood for the Cartesian label a (i.e. v ia u a = 
>~2 a Vi a u a ) and where i>b{n) is the free energy in the bulk. We conveniently choose 
the double well form 

Mn) = Pc (f„ + I) 2 {vl - 2z/ n + 3 - 2/3r w ) (2) 

where v n = (n-n c )/n c , t w = (T c -T)/T c mdp c = 1/8, n c = 7/2 and T c = 4/7 
are the critical pressure, density and temperature respectively and (3 is a constant 
typically equal to 0.1. 

The derivative term in equation ([TJ models the free energy associated with 
density gradients at an interface, k is related to the surface tension 7 by Q 
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7 = -^p~ c {(3T w f' 2 n c . (3) 

^c( n s) = —4>i n s is the Cahn surface free energy [ 8 1 which controls the wetting 
properties of the fluid. In particular (pi can be used to tune the contact angle. 

2.2 Navier-Stokes equations 

The dynamics of the droplet is described by the Navier-Stokes equations for a 
non-ideal gas 

d t (nu a ) + d[3(nu a up) = -d p P a p + vdp [n(dpu a + d a up + 5 a/ 3<9 7 u 7 )] + nF a , 
d t n + d a (nu a ) = (4) 

where u(r) is the fluid velocity, v the kinematic viscosity and F a gravitational 
field. The pressure tensor P a p is related to the free energy by [9] 

dT 

d{d a n) 

where T = ipb — ^bn + K,(d a n) 2 /2 and fi b = Ap c (l — (3t w ) / n c is the bulk chemical 
potential. 
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Figure 1: Topology of a D3Q15 lattice. The directions i are numbered and corre- 
spond to the velocity vectors v$. 

2.3 The lattice Boltzmann algorithm 

We solve the equations of motion © by using a lattice Boltzmann algorithm. This 
approach follows the evolution of partial distribution functions fa on a regular, d- 
dimensional lattice formed of sites r. The label i denotes velocity directions and 
runs between and z. DdQz + 1 is a standard lattice topology classification. 
The D3Q15 lattice we use here has the following velocity vectors Vj! (0, 0, 0), 
(±1, ±1, ±1), (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) in lattice units as shown in fig. □ 
The lattice Boltzmann dynamics are given by 

fi(r + Ahn, t + At) = fi(r, t) + - (jf (r, t) - /«(r, £)) + n W<T v ia F a (6) 

T 

where At is the time step of the simulation, r the relaxation time, a labels ve- 
locities of different magnitude, w\ = 1/3, W2 = 1/24. f^ 9 is the equilibrium 
distribution function which is a function of the density n = J2i=o fi an( ^ me ^ ul< ^ 
velocity u, defined through the relation 



nu = 

i=0 



The relaxation time tunes the kinematic viscosity as 



f Ar) 2 1 1 
u = i=^-i(r - - (8) 

At 3 V 2 J ' J 



where Ar is the lattice spacing. 
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It can be shown [ 10 1 that equation © reproduces the Navier-Stokes equations 
of a non-ideal gas @ if the local equilibrium functions are chosen as 



f- q = A a + B a u a v ia + C a u 2 + D a u a Ui3V ia Vii3 + G aa pv ia Vip, i > 0, 

fo g = (9) 



i=l 

A possible choice of the coefficients is [fTTI 



A<T = (p° ~ ^(^ n ) 2 ~ Knd aa n + vu a d a n 
R w CT n w a n 3w rr n 



C 



2c 2 ' 2c 4 



^177 = (k(<9 7 7i) 2 + 2i>u 1 d 1 n) , G277 = 0, 
G 2 ~/6 = t^-t {K(d~n)(d s n) + v(u 7 d s n + u s d 7 n)) (10) 
where c = Ar/At and p = - if>b = Pc(v n + 1) 2 (3^ -2u n + l- 2(3t w ). 

2.4 Wetting boundary conditions 

The major challenge in dealing with patterned substrates is to handle the boundary 
conditions correctly. We consider first wetting boundary conditions which control 
the value of the density derivative and hence the contact angle. For flat substrates 
a boundary condition can be established by minimising the free energy CQ> 1 8 1 

s-Vn = -— (11) 

K 

where s is the unit vector normal to the substrate. It is possible to obtain an 
expression relating (pi to the contact angle 9 as 



0i = 2f3r w y/2p c K sign ~ °) y cos f ( X ~ cos f ) ( lT > 

where a = cos _1 (sin 2 6) and the function sign returns the sign of its argument. 

Equation ([HI) is used to constrain the density derivative for sites on a flat part 
of the substrate. However, no such exact results are available for sites at edges or 
corners. We work on the principle that the wetting angle at such sites should be 
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Figure 2: Sketch of a post on a substrate. Encircled numbers label sites in different 
topological positions. Labels 26 and 27 denote sites on the bottom (z = z min ) and 
the top (z = z max ) of the domain respectively. 

constrained as little as possible so that, in the limit of an increasingly fine mesh, it 
is determined by the contact angle of the neighbouring flat surfaces. 

For edges (labels 9 — 12 in fig. O and corners (labels 1 — 4) at the top of the 
post each site has 6 neighbours on the computational mesh. Therefore these sites 
can be treated as bulk sites. 

At bottom edges where the post abuts the surface (labels 13 — 16 in fig. EJ) 
density derivatives in the two directions normal to the surface (e.g. x and z for 
sites labeled 13) are calculated using 



where the middle term constrains the density derivative in the appropriate direc- 
tion xory. 

At bottom corners where the post joins the surface (labels 5 — 8 in fig. |2) 
density derivatives in both the x and y directions are known. Therefore these sites 
are treated as planar sites. 




(13) 
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2.5 Velocity boundary conditions 

We impose a no-slip boundary condition on the velocity by determining the miss- 
ing fields which fulfill the no- slip condition given by equation © with u = 0. 
This does not uniquely determine the fa's. For most of the cases (i.e. 1 — 20) 
arbitrary choices guided by symmetry are used to close the system. This is no 
longer possible for sites 21 — 27 where four asymmetrical choices are available. 
Selecting one of those solutions or using a simple algorithm which chooses one of 
them at random each time step leads to very comparable and symmetrical results. 
Hence we argue that an asymmetrical choice can be used. Possible conditions, 
which are used in the results reported here, are listed in appendix A. 

The conservation of mass is ensured by setting a suitable rest field, /o, equal 
to the difference between the density of the missing fields and the one of the fields 
entering the solid after collision. 

In a hydrodynamic description of wetting contact line slip must be introduced 
in some way. As with other phase field models slip appears naturally within the 
lattice Boltzmann framework. The mechanism responsible for the slip is evapora- 
tion and condensation of the fluid because of a non-equilibrium curvature of the 
contact line 11710. 

3 Results 

We consider the superhydrophobic behaviour of droplet spreading on a substrate 
patterned by square posts arranged as in fig. |3] The size of the domain is L x xL y x 
L z = 80 x 80 x 80 and the height, spacing and width of posts are h = 5, d = 8 
and w = 4 respectively. A spherical droplet of radius R = 30 is initially centered 
within the domain and just touches the post tops. The contact angle 6 eq = 110° 
is set on every substrate site. The surface tension and the viscosity are tuned by 
choosing parameters k = 0.002 and r = 0.8 respectively. The liquid density ni 
and gas density n g are set to n\ = 4.128 and n g = 2.913 and the temperature 
T = 0.4. 

Simulation and physical parameters can be related by choosing a length scale 
Ar, a time scale At and a mass scale Am. A simulation parameter with di- 
mensions [L] ni [T] n2 [M] n3 is multiplied by Ar ni At n2 Am™ 3 to give the physical 
value. For example, considering a 1 mm droplet of a fluid of kinematic viscosity 
v — 3 • 10~ 5 m 2 s _1 and surface tension 1 • 10~ 3 Nnr 1 leads to Ar = 1.7 • 10 -5 
m, At = 9.3 ■ 10~ 7 s, Am = 1.6 • 10~ 12 kg. That implies n x = 1.4 • 10 3 kgm" 3 . 
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Figure 3: Plan view of the substrate. Shaded areas are posts. 

However, as with all diffuse interface models the liquid-gas density difference is 
unphysically small and the width of the interface is too large. This must be taken 
into account by renormalising the time scale. 

3.1 Equilibrium states 

Fig. HJshows the final state attained by the droplet for different substrates or initial 
conditions. For comparison fig. Uta) shows a planar substrate. The equilibrium 
contact angle is 9^ lat = 110° = 6 eq as expected ifTTl . In fig. |4[b) the substrate 
is patterned and the initial velocity of the drop is zero. Now the contact angle is 
9 s = 156°, a demonstration of superhydrophobic behaviour. Fig. Etc) reports an 
identical geometry but a drop with an initial impact velocity. Now the drop is able 
to collapse onto the substrate and the final angle is 9 C = 130°. These angles are 
compatible with the ones reported in Q. 

Superhydrophobicity occurs over a wide range of d, the distance between the 
posts. For suspended drops of R = 30 and d> 12 the drop resides on a single post 
and the contact angle is 170°. For d < 12 the contact angle lies between 148° and 
156° with the range primarily due to the commensurability between drop radius 
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Figure 4: Final states of a spreading droplet. The right column reports vertical 
cuts across the centre of the drop, (a) The substrate is flat and homogeneous, (b) 
The substrate is decorated with posts and the initial velocity of the droplet is zero, 
(c) Same geometry as (b) but the droplet reaches the substrate with a velocity 
u z = -0.01. 
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and post spacing. 



3.2 Kinetics of the suspended to collapsed transition 

We now investigate the kinetics of the transition between the suspended and col- 
lapsed droplet states. For the parameter values used in the simulations presented 
in fig. |U the state with the drop suspended on the posts has a slightly higher free 
energy than the collapsed state. However as the drop penetrates the grooves the 
area of the contact between liquid and solid increases. Because the substrate is 
hydrophobic this creates a free energy barrier hindering the transition from the 
suspended to the collapsed states. Work must be provided, by an impact velocity 
or gravity say, to allow the transition to proceed lfT3l l3ll. 

We follow the transition pathway by considering a spherical droplet of radius 
R = 30 initially centered within the domain and just touching the top of the 
posts. A gravitational field F z is turned on at time t = 70000, and turned off at 
t = 200000. Fig. 13 shows cross sections of the drop as it undergoes the transition 
from the suspended to the collapsed state for F z = — 5 • 10~ 7 . Note how the 
substrate interstices are filled successively from the drop centre to its edges. 

The time evolution of the total free energy of the system is presented in fig.® 
for two different values of F z . We plot \P T the total free energy; the volume ty v and 
surface contributions to the free energy which arise from the first and second 
integral in (1) respectively and the contribution from the gravitational force 



The solid lines in fig. ® correspond to the snapshots in fig. [3 After gravity is 
switched on the drop is pushed down and hence ty g decreases. On the other hand 
\l/ s increases as the surface is covered by liquid rather than by gas corresponding 
to a free energy gain for a hydrophobic surface. ty v also initially increases as the 
interface is deformed. Once the drop touches the bottom of the interstices 
drop sharply because parts of the interface have just vanished. At the same time 
\l/ s grows significantly because a part of the surface is now in contact with the 
liquid rather than with air. At t = 200000, the gravitational field is turned off. 
The squeezed droplet dewets to recover a spherical shape and ^ s decreases. The 
total free energy is slightly smaller in the collapsed state. 

The figures show that both the surface and the volume energies increase during 
the transition. Hence, as the total free energy must decrease, the transition could 




(14) 



10 



70000 



70r 
60 
50- 
40- 




n n n R^^^n n n 



75000 



70 p 

60 

50- 



In n n^rmrrmti n nl 



100000 



n n-n-n-n-n-n" 



U = 124000 



129500 




138000 




| | j 1 ^ | | | | | | | | | | | | ^pnj | | 
20 40 60 80 



140000 




| | p~| | | | | | | | | | | | | | | 



300000 



| | p~| | || || || || || | |~~] | | 



Figure 5: Transition from a suspended to a collapsed state. The gravity field 
F z = 5- 10~ 7 . These cuts are vertical cross sections across the centre of the domain 
where the dark gray areas are the posts and the pale ones are liquid regions. 
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Figure 6: Evolution of the free energies of superhydrophobic drops. A body force 
F z is turned on at t = 70000 and turned off at t = 200000. Solid and dashed lines 
corresponds to free energies when F z = — 5- 10~ 7 and F z = — 3- 10~ 7 respectively. 
The axes are in simulation units with free energy chosen to be zero at equilibrium. 
In each case the larger diagram highlights the free energy and time regimes where 
the most interesting behaviour occurs. The insets show the behaviour throughout 
the simulation. 
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not occur without the addition of a gravitational force. The decrease in the gravi- 
tational free energy offsets the increase in the surface and volume terms allowing 
the transition to proceed. Essentially the addition of F z removes the barrier which 
the total free energy must surmount. 

A lower gravitational field, F z = 3 ■ 10~ 7 , does not allow the drop to overcome 
the free energy barrier and it remains in the suspended state as shown by the 
dashed lines in fig. 

For example, considering again a 1 mm droplet of a fluid of kinematic viscos- 
ity v = 3 ■ 10~ 5 m 2 s _1 and surface tension 1 • 1CT 3 Nm _1 , F z = 5 • 1CT 7 would 
correspond to 9.7 ms~ 2 in physical units. 

4 Conclusion 

We have used a lattice Boltzmann approach to solve the equations of motion de- 
scribing the dynamics of a drop on topologically patterned substrate. The ap- 
proach allows us to simulate the dynamic and equilibrium properties of a drop 
with a given size, surface tension, contact angle and viscosity. Because interfaces 
appear naturally within the model it is particularly well suited to looking at the 
behaviour of evolving drops. 

In particular we have considered a droplet positioned on an array of posts and 
shown that it is possible to reproduce the superhydrophobic behaviour seen in 
experiments. The ability of the algorithm to follow the drop kinetics has enabled 
us to investigate the transitions between the suspended state where the drop lies 
on top of the posts and the collapsed state where it fills the spaces between them. 
We find that the substrate interstices are filled successively starting from the drop 
centre. 

There are many avenues open for further investigation. For example we are 
currently investigating how drops move across hydrophobic surfaces to compare 
to recent experiments. It would also be interesting to follow the spreading of drops 
on surfaces with topological imperfections, a problem of concern in the quality of 
images formed in ink-jet printing. 
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A Possible boundary conditions 

We list here the boundary conditions used to define missing distribution functions, 
i.e. those that stream from positions outside the simulation box. 



Label (see fig. 2) 


Conditions 


1 


fl3 = fu 


2 


h — h 


3 


h — ho 


4 


fll = fl2 


5 


h = fe 

fis = (h ~ h ~ h + h)/2 + h + fu - fio 

fll = (fl - h)/2 -fg + f W + fl2 
f7 =(-/3 + /4)/2 + / 8 -/ 9 + /l0 


6 


/l3 = (/ 3 - /4)/2 - fll + fl2 + fu 

h = (fi - h)/2 - f U + f 10 + f l2 

h = (-h + h-h + h)/2 + f 8 + fii - fu 


7 


fn = (h - h)/2 - f 13 + f 12 + fu 

h = (fi -h-h + h)/2 + As - fu + ho 

fl =(-/l+/ 2 )/2 + / 8 -/i 3 + /l4 


8 


hi =(h-h + h-h)/2 + h-h + h 2 
h = (-h + / 4 )/2 - h + h + ho 
fn = (-h + h)/2 -h + h + fu 


9 


fi3 — fu ; h = fs 


10 


h =ho ; h = h 


11 


fo = flO fll = fl2 


12 


fi3 = fu ; fu = fu 


13 


y 5 = f & 

h = 2(-f w + f 9 + f u - f 12 ) + h 
fi3 = (h ~ h)/2 - fu + fl2 + fu 

h =(-h + h)/2 + h-h + ho 
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J (UL/wl 1 jvv 11. . . 1 


Conditions 


1 A 

14 


is 


— 


fa 




fa 


— 


(fa ~ fa)/ 2 - fu + /iQ + /l2 




fa 


— 


(-fa + / 2 )/2 + fa - fa 3 + fu 




fa 


= 


2(-/l2 + fu + fl3 ~ fu) + fa 


15 


fa 


— 


fa 




fa 


— — 


2(-/i4 + A + /is - fa) + A 




£ 

fll 


— 


[fa - fa)/ 2 - fl3 + fl2 + fu 




£ 

fa 


— 


( £ i £ \ 1 C\ i £ £ i £ 

{-fa + /4J/2 + fa- fa + fio 


lo 


fa 


— — 


fa 




fu 


— — 


(fa - fa)/ 2 ~ fa + fio + fl2 




£ 

J13 




( £ 1 £ \ /O £ 1 £ 1 £ 

(-fa + /2J/2 - fa + fa + fu 




fa 


— 


n / £ 1 £ 1 £ £ \ 1 £ 

2(-/io + fa + fa- fa) + fa 


1 7 
1 I 


£ 

ho 


— 


£ £ £ 

j9 ; J13 — J14 


lo 


fa 


— 


-f £ £ 

J 8 ; J12 = jn 




fa 


= 


/10 ; /l4 = /l3 


on 
2U 


fa 


— 


A ; /n = /i2 


21 


fa 


— 


A ; A — A 




fl2 


— 


(-/a + A)/2 + fu 




fl3 


— 


(~fa + A)/2 + /14 




fio 


— 


(fa - A + A - A)/2 + A 


22 


fa 


— 


A ; A = A 




fa 


— 


(-fa + A)/2 + Ao 




fu 


— 


(fa-fa + fa- /e)/2 + /i 3 




fl2 


— 


(-fa + A)/2 + Ai 


23 


fa 




A ; A = A 




fll 




(A-AV2 + A2 




fu 




(A-A)/2 + A 3 




fa 




Ao + (-A + A - A + A)/2 


24 


fa 




A ; A — A 




fll 




(A-A)/2 + A 2 




fio 




(A-A)/2 + A 




fl3 




(-A + fa - fa + fa)/2 + fa A 


25 


fa 




fa ', fa — fa 




fll 




(fi-fa + fa-fa)/2 + fa 2 
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Label (see fig. 2) 


Conditions 




h =(-/ 3 + /4)/2 + /io 

Jl3 - {-Jl + /2J/2 + /14 


26 


A — fs ; /s — /6 

/ll =(/l-/ 2 + /3-/4)/2 + /i2 

/ 9 =(-/ 3 + /4)/2 + /io 

/l3 =(-/l + / 2 )/2 + /i4 


27 


/6 = /s ; /s = A 
/10 =(/ 3 -/4)/2 + / 9 

/l4 =(/l-/ 2 )/2 + /l 3 

/l2 =(-/3 + /4-/l + / 2 )/2 + /n 
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